{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Handling uncertainty with quantile regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-04T17:52:47.886294Z",
     "iopub.status.busy": "2023-12-04T17:52:47.885937Z",
     "iopub.status.idle": "2023-12-04T17:52:48.106622Z",
     "shell.execute_reply": "2023-12-04T17:52:48.106311Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Quantile regression](https://www.wikiwand.com/en/Quantile_regression) is useful when you're not so much interested in the accuracy of your model, but rather you want your model to be good at ranking observations correctly. The typical way to perform quantile regression is to use a special loss function, namely the quantile loss. The quantile loss takes a parameter, $\\alpha$ (alpha), which indicates which quantile the model should be targeting. In the case of $\\alpha = 0.5$, then this is equivalent to asking the model to predict the median value of the target, and not the most likely value which would be the mean. \n",
    "\n",
    "A nice thing we can do with quantile regression is to produce a prediction interval for each prediction. Indeed, if we predict the lower and upper quantiles of the target then we will be able to obtain a \"trust region\" in between which the true value is likely to belong. Of course, the likeliness will depend on the chosen quantiles. For a slightly more detailed explanation see [this](https://medium.com/the-artificial-impostor/quantile-regression-part-1-e25bdd8d9d43) blog post.\n",
    "\n",
    "As an example, let us take the [simple nowcasting model we built in another notebook](building-a-simple-nowcasting-model.md). Instead of predicting the mean value of the target distribution, we will predict the 5th, 50th, 95th quantiles. This will require training three separate models, so we will encapsulate the model building logic in a function called `make_model`. We also have to slightly adapt the training loop, but not by much. Finally, we will draw the prediction interval along with the predictions from for 50th quantile (i.e. the median) and the true values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-04T17:52:48.108477Z",
     "iopub.status.busy": "2023-12-04T17:52:48.108356Z",
     "iopub.status.idle": "2023-12-04T17:52:48.657556Z",
     "shell.execute_reply": "2023-12-04T17:52:48.656135Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy0AAAIQCAYAAACMg4HBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABe/klEQVR4nO3deXxU1f3/8fedCdk39kUCqCjiglVBiq0LigKiRa3FIgqoP9qquNRqldZWsVqtdaPfttDvVwVbQEGt4JeqfCkFBUTEBcUNkUIBZd8SAlnm3vP7I8ww+0ySuZMweT0fjyhz587cc+cQzTvnfM6xjDFGAAAAANBMeZq6AQAAAAAQD6EFAAAAQLNGaAEAAADQrBFaAAAAADRrhBYAAAAAzRqhBQAAAECzRmgBAAAA0KwRWgAAAAA0a4QWAAAAAM0aoQUAAABAs0ZoAYAg06ZNk2VZsixLS5cujXjeGKOysjJZlqVLLrkk6nvs3btXubm5sixLn3/+edRzxo4dG7hO+Fdubm6D2v7NN9/ommuuUa9evVRUVKTS0lKdeeaZeu6552SMCTm3R48eMa9/3HHHxbzG0qVLA+ft3LkzqXY5jqNHH31URx99tHJzc9WnTx89//zzcV9TW1urE088UZZl6bHHHmvwfUazcuVKjR8/XieddJIKCgrUrVs3jRgxQl9++WXM9k+ePFnf+ta3lJeXp7Zt2+r888/XRx99FPMaM2bMkGVZKiwsjHivadOm6Xvf+57KyspUUFCgk08+WQ8++KCqqqoi3mfbtm267rrr1KFDB+Xl5en000/Xiy++mPAeL7zwQlmWpfHjxyc8FwCOBFlN3QAAaI5yc3M1c+ZMffe73w05/uabb2rz5s3KycmJ+doXX3xRlmWpU6dOmjFjhh588MGo5+Xk5Ojpp5+OOO71ehvU5p07d2rz5s268sor1a1bN9XW1mrBggUaO3as1qxZo9/+9reBc5966int378/5PX/+c9/dO+99+qiiy6K+v6O4+iWW25RQUGBKisrk27XL3/5Sz3yyCMaN26c+vXrp7lz5+rqq6+WZVn64Q9/GPU1//Vf/6WNGzc2+j6j+d3vfqdly5bpBz/4gfr06aOtW7fqj3/8o04//XS98847Ovnkk0POv/766zVjxgyNHj1a48ePV2VlpT788ENt37496vvv379fP//5z1VQUBDx3IEDB3Tdddfp29/+tn7yk5+oQ4cOWr58ue677z4tXLhQ//rXv2RZliSpvLxc3/3ud7Vt2zbddttt6tSpk2bPnq0RI0ZoxowZuvrqq6Ne/+9//7uWL18e9zMAgCOOAQAETJ061UgyV1xxhWnXrp2pra0NeX7cuHHmjDPOMN27dzfDhg2L+h7nnHOOueKKK8xPf/pTc/TRR0c9Z8yYMaagoCDl7Y/mkksuMQUFBcbn88U97ze/+Y2RZJYtWxb1+cmTJ5u2bdua2267zUgyO3bsSHjtzZs3m1atWpmbb745cMxxHHP22Webrl27Rm3Ttm3bTElJiXnggQeMJPP73/8+4XWMSf4+ly1bZqqrq0OOffnllyYnJ8eMGjUq5PisWbOMJPP3v/89qTYYY8zdd99tevXqZUaNGhXRx9XV1VE/34kTJxpJZsGCBYFjjz76qJFkFi5cGDhm27bp16+f6dSpU8Q9GGPMwYMHTY8ePQKfXfDnDgBHMqaHAUAUI0eO1K5du7RgwYLAsZqaGr300ksxf8MtSRs3btSSJUv0wx/+UD/84Q+1fv16vf32241qy7p167Ru3boGv75Hjx46cOCAampq4p43c+ZMHX300TrrrLMintu9e7fuvfdePfDAAyotLU362nPnzlVtba1uuummwDHLsnTjjTdq8+bNUUcE7rnnHvXq1UvXXHNN0teRkr/Ps846S9nZ2SHHjjvuOJ100kkR0/meeOIJnXnmmbr88svlOE7CEaa1a9fqySef1BNPPKGsrMjJDNnZ2VE/38svv1ySQq6/ZMkStW/fXueff37gmMfj0YgRI7R161a9+eabEe/z6KOPynEc3XnnnXHbCQBHGkILAETRo0cPDRgwIKT24vXXX9e+fftiTmmSpOeff14FBQW65JJLdOaZZ+rYY4/VjBkzYp6/c+fOiK/y8vKQcy644AJdcMEFSbf94MGD2rlzpzZs2KDnnntOU6dO1YABA5SXlxfzNR9++KE+//zzmIHsV7/6lTp16qQf//jHSbfD/74FBQXq3bt3yPEzzzwz8Hywd999V88995yeeuqpwDSpWBpyn7EYY7Rt2za1a9cucKy8vFzvvvuu+vXrp1/84hcqKSlRYWGhjjnmGM2ePTvq+9x+++0aOHCgLr744npdf+vWrZIUcv3q6uqo95Kfny9Jev/990OOb9y4UY888oh+97vfNegzAIDmjNACADFcffXVmjNnjg4ePCiprrj63HPPVZcuXWK+ZsaMGRo+fHjgh8arrrpKs2fPls/nizi3srJS7du3j/gaMWJEo9o9adIktW/fXkcffbTGjh2rb3/723rhhRfivsYfrEaNGhXx3Mcff6y//OUveuKJJ+pdb7NlyxZ17NgxIoB07txZUl1RvZ8xRrfccouuuuoqDRgwIOF7N+Q+Y5kxY4a+/vprXXXVVYFj69atkzFGL7zwgp599lk9+uijmjFjhtq3b68f/vCHeuONN0Le4x//+If+7//+T0888US9r//oo4+quLhYQ4cODRzr1auXNm/erP/85z8h5y5ZskSS9PXXX4cc/9nPfqbTTjstbqgGgCMVhfgAEMOIESN0++23a968eRoyZIjmzZunP/zhDzHP//jjj7V69Wo9/PDDgWMjR47Ub3/7W82fP1/Dhg0LOT83N1f/+7//G/E+wb9tl6QNGzbUq90jR45U3759tWPHDs2bN0/btm0LBK9oHMfRCy+8oNNOOy1iRESSbr31Vg0dOjRmgX48Bw8ejLpogX+FtOB2TZs2TatXr9ZLL72U1HvX9z5j+eKLL3TzzTdrwIABGjNmTOC4f6GCXbt26Z133lH//v0lSd/73vd09NFH68EHH9SQIUMk1U0d/OlPf6qf/OQnOvHEE+t1/d/+9rf65z//qT//+c8hU+/+3//7f5oyZYpGjBihJ598Uh07dtTs2bP1yiuvSAr97BYtWqSXX35ZK1asqPf9A8CRoMlCy1tvvaXf//73ev/997Vlyxa98soruuyyy5J+/f3336+JEydGHM/Pz6/XqjYAEEv79u01aNAgzZw5UwcOHJBt27ryyitjnj99+nQVFBTomGOO0VdffSWp7ofzHj16aMaMGRGhxev1atCgQSlvd/fu3dW9e3dJdT/Y/+hHP9KgQYO0Zs2aqNOG3nzzTX399df66U9/GvHcrFmz9Pbbb+uTTz5pUFvy8vJUXV0dcdy/vK+/PeXl5ZowYYLuuusulZWVJfXe9b3PaLZu3aphw4appKREL730UshIkv89jj766EBgkaTCwkJdeumlmj59unw+n7KysvTkk09q586dUf+/FM+sWbN077336oYbbtCNN94Y8lyfPn00c+ZM/eQnP9F3vvMdSVKnTp301FNP6cYbbwwsp+zz+XTrrbfq2muvVb9+/ep1fQA4UjTZ9LDKykqdeuqp+tOf/tSg1995553asmVLyNeJJ56oH/zgByluKYCW7Oqrr9brr7+uKVOmaOjQoTGL0I0xev7551VZWakTTzxRxx13XOBrw4YNmjt3bsQSw+ly5ZVXatOmTXrrrbeiPj9jxgx5PB6NHDky4rm77rpLP/jBD5Sdna0NGzZow4YN2rt3ryRp06ZNIdO7ouncubO2bt0asX/Kli1bJCkw1e6xxx5TTU2NrrrqqsB1Nm/eLEnas2ePNmzYkLDAPtF9htu3b5+GDh2qvXv36o033oiY9ud/3LFjx4jXdujQQbW1taqsrNS+ffv04IMPaty4cSovLw+0f//+/TLGaMOGDVGXR16wYIFGjx6tYcOGacqUKTHv6ZtvvtG7776r5cuX6z//+Y+OOeYYSdLxxx8vSfrrX/+qNWvW6Mc//nHg2v7RuYqKCm3YsEEHDhxI6jMBgGarKZcu85NkXnnllZBjVVVV5mc/+5np0qWLyc/PN2eeeaZZtGhRzPdYtWqVkWTeeustdxsLIKP5lzxeuXKlMcaYiooKk5eXZySZWbNmBc4LX/J40aJFRpJ54IEHzIsvvhjy9d///d9Gkvnb3/4WOD+dSx7PmTMnov1+VVVVprS01Jx//vlRXysp7tepp54a99p//OMfjSTz6aefhhyfMWNGyH+zx4wZk/BaH374YYPvM9zBgwfN2WefbfLz883bb78d87xOnTqZsrKyiOPXXnutyc3NNbZtm/Xr1yds+/Dhw0Ne/84775iCggJz1llnmQMHDiRsb7C77rrLSDJr1qwxxhhz3333Jbx++P9jAeBI02xrWsaPH6/PPvtML7zwgrp06aJXXnlFQ4YM0erVq6Pu1vz000/r+OOP19lnn90ErQWQqQoLCzV58mRt2LBBl156aczz/FPD7rrrrqg72v/+97/XjBkz6r2Mr6TAcsfHHnts3PN27Nih9u3bRxx/5plnZFmWTj/99IjnXnvtNe3duzdqAb6kQP1EsBdeeEGzZs3SX//6V3Xt2jVum4YPH66f/vSn+vOf/6w//vGPkupGpaZMmaKjjjoqsPzvrbfeGjFFePv27frxj3+ssWPHavjw4Tr66KPrfZ/+Fdm6desWWHXLtm1dddVVWr58uebOnRu36P+qq67SpEmTtGDBAl144YWB95w7d67OP/98eTwedejQIern9Ic//EHLly/X888/H1h4QKpb1njYsGHq0aOH5s2bV6+VvtauXaspU6bokksuCYy0/PCHP9S3vvWtiHMvv/xyXXzxxRo3blzI9DYAOBI1y9CyceNGTZ06VRs3bgwMz99555164403NHXq1IjdjquqqjRjxgzdc889TdFcABkuuDg7murqar388su68MILowYWqa54e9KkSdq+fbs6dOggqa4WYfr06VHPv/zyywM7qvuXO05UkP/QQw9p2bJlGjJkiLp166bdu3fr5Zdf1sqVK3XLLbeoZ8+eEa+ZMWOGcnJy9P3vfz/qe0arNVy1apUkaejQoSGLBixevFgDBw7Ufffdp/vvv1+S1LVrV91+++36/e9/r9raWvXr109z5szRkiVLNGPGjEANyemnnx4Rqvz3e9JJJ4W0oz73+cc//lETJ07UokWLdN5550mqW2Xr1Vdf1aWXXqrdu3dH9EFwsJwwYYJmz56t73//+7rjjjtUUlKiKVOmqLa2NvD/ovz8/Kif05w5c/Tuu++GPFdRUaHBgwdrz549uuuuu/SPf/wj5DXHHntsSIjyT3vu1q2b1q9fr8mTJ6tNmzYh08lOOOEEnXDCCRHXl+rqcepTLwoAzVWzDC2rV6+WbduB3yL5VVdXq23bthHnv/LKK6qoqEj4gwUAuOEf//iH9u7dG3ck5tJLL9Xjjz+uF154Qbfeequkuv+mXXvttVHPX79+fSC0JGvYsGFat26dnn32We3YsUO5ubnq06ePpk6dGvW/j+Xl5frHP/4RKERvLH/NTvCogiQ98sgjat26tf7yl79o2rRpOu644zR9+vS4m3TGU9/7DOcPXf/7v/8bdfW24NDSsWNHLV26VHfeeaeefPJJ1dbWasCAAZo+fbpOPfXUerd9165d2rRpkyRF/UXbmDFjQkLLqaeeqqlTpwb2kBkxYoQmTpwYCL4A0FJYxoRVRzZFIywrZPWwWbNmadSoUfr0008j9gQoLCxUp06dQo5dcMEFKi4ujjo8DwBIj5///Od6/vnn9dVXX0Vd5hgAgIZqliMtp512mmzb1vbt2xPWqKxfv16LFi3Sq6++mqbWAQCiWbRokX71q18RWAAAKddkoWX//v2BfQykuvCxatUqtWnTRscff7xGjRql0aNH6/HHH9dpp52mHTt2aOHCherTp0/IXgfPPvusOnfuHLKLMAAg/VauXNnUTQAAZKgmmx7mL9gMN2bMGE2bNk21tbV68MEH9de//lVff/212rVrp29/+9uaOHGiTjnlFEl1uzh3795do0eP1kMPPZTuWwAAAACQBs2ipgUAAAAAYvE0dQMAAAAAIB5CCwAAAIBmLe2F+I7j6JtvvlFRUZEsy0r35QEAAAA0E8YYVVRUqEuXLvJ4Yo+npD20fPPNNyorK0v3ZQEAAAA0U5s2bVLXrl1jPp/20FJUVCSprmHFxcXpvnyL5fP59N5776lv377KymqW2/OgAejXzES/Zi76NjPRr5mJfk2P8vJylZWVBTJCLGnvAf+UsOLiYkJLGvl8PhUUFKi4uJhvvAxCv2Ym+jVz0beZiX7NTPRreiUqG6EQHwAAAECzRmgBAAAA0KwRWgAAAAA0a0zQAwAAOELYtq3a2tqmbkaL4PP5ZFmWqqqqqGlphFatWsnr9Tb6fegBAACAZs4Yo61bt2rv3r1N3ZQWwxij/Px8bdy4kb0FG6m0tFSdOnVq1OdIaAEAAGjm/IGlQ4cOys/P54foNDDG6MCBA3zejeD/DLdv3y5J6ty5c4Pfi9ACAADQjNm2HQgsbdu2bermtBjGGNm2rdzcXEJLI+Tl5UmStm/frg4dOjR4qhiF+AAAAM2Yv4YlPz+/iVsCNIz/725j6rEILQAAAEcAftuPI1Uq/u4SWgAAAAA0a4QWAAAAHLEWL14sy7JYWS3DEVoAAACQUpZlxf26//77G/S+5513nm6//faUthVHBlYPAwAAQEpt2bIl8OdZs2bp17/+tdasWRM4VlhYGPizf5UuNnBEPPUaabn//vsjkvIJJ5zgVtsAAAAQxjGOyn0HmuzLMU7CNnbq1CnwVVJSIsuyAo+/+OILFRUV6fXXX9cZZ5yhnJwcLV26VGPHjtVll10W8j633367zjvvPEnS2LFj9eabb2rSpEmBn0M3bNgQOPf9999X3759lZ+fr7POOiskJOHIV+9Ie9JJJ+mf//zn4TcgFQMAAKTNfrtKP147pcmu/5fjfqLirMYvv3zPPffoscce0zHHHKPWrVsnPH/SpEn68ssvdfLJJ+uBBx6QJLVv3z4QXH75y1/q8ccfV/v27fWTn/xE119/vZYtW9bodqJ5qHfiyMrKUqdOndxoCwAAAFqIBx54QBdeeGHS55eUlCg7O1v5+flRfxZ96KGHdO6550qqC0TDhg1TVVWVcnNzU9ZmNJ16F+KvXbtWXbp00THHHKNRo0Zp48aNbrQLAAAAGaxv374pfb8+ffoE/ty5c2dJdbuwIzPUa6Slf//+mjZtmnr16qUtW7Zo4sSJOvvss/XJJ5+oqKgo6muqq6tVXV0deFxeXi5J8vl88vl8jWg66sO27UChGzIH/ZqZ6NfMRd9mJrf71efzyRgT8tWU6tsG/7nh/87Pzw95H8uyIt67pqYm5DXRru//c1ZWVkS7/H3TEM3l884E/s8x2s//yeaBeoWWoUOHBv7cp08f9e/fX927d9fs2bN1ww03RH3Nww8/rIkTJ0Ycf++991RQUFCfyzcZ5+BBWV6vrOzspm5KgzmOo4qKCr377rvyeFjpOlPQr5mJfs1c9G1mcrtfLctSfn6+Dhw4INu2ZRlHj3e5NuXXSbo9VbYqrcqkz/f/8rqysu41VVVVgcetWrUKnFdaWqrVq1cHzpOkDz74QK1atQoc83q9qq6uDjkn2vsdPHhQknTgwIGQc+vDGCPHcVRZWZmSHd1bsurqatXU1Ojjjz+OCIHJ9k+jquhLS0t1/PHH66uvvop5zoQJE3THHXcEHpeXl6usrEx9+/ZVcXFxYy6fNjUffyBPUbGyju7Z1E1pMNu2tXLlSvXr109er7epm4MUoV8zE/2auejbzOR2v1ZVVWnjxo3Kz88P1GcUKfoMl+YoJydHkgK/rPbfQ0FBQcgvsAcPHqxJkybp5Zdf1oABAzR9+nR9/vnnOu200wLnHXPMMfrggw+0Y8cOFRYWqk2bNlHfLy8vT1LdaE5Df0lujFFlZaUKCgoILY3k9XqVnZ2tnj17RtQY+WdhJdKo0LJ//36tW7dO114bO+3n5OQE/rKGXDgr64hZeczn2PI4R/764ZZlyev1HvH3gVD0a2aiXzMXfZuZ3OzXrKyskO0mjjT+Nkf7d/D9DBkyRL/61a909913q6qqStdff71Gjx6t1atXB8676667NGbMGJ100kk6ePCg1q9fH/X9Yl2jIW0/Uj/35sT/GUb7+T/Z75l6fWfdeeeduvTSS9W9e3d98803uu++++T1ejVy5Mj6vM2RxxiZoLocAAAAJGfs2LEaO3Zs4PF5550Xs05k4sSJUcsK/I4//ngtX7485FiPHj0i3u9b3/oWtSgZpl6hZfPmzRo5cqR27dql9u3b67vf/a7eeecdtW/f3q32NQ+OkamuaupWAAAAAC1SvULLCy+84FY7mjfjyFSz0hkAAADQFFi6JEmmppphRgAAAKAJEFqSYYxk2xL7ygAAAABpR2hJgjFGxrZlfLVN3RQAAACgxSG0JIORFgAAAKDJEFqSYRzJtmUILQAAAEDaEVqS4RjJsSWmhwEAAABpR2hJhjEyto+RFgAAAKAJEFqSQk0LAABAczV27FhddtllgcfnnXeebr/99ka958CBA3X33Xc3rmFIGUJLMoyRsR1WDwMAAKiHsWPHyrIsWZal7Oxs9ezZUw888IB8Lv8i+O9//7t+85vfJHXu4sWLZVmW9u7dG3L85Zdf1r333utC69AQWU3dgCOBMZJsHyMtAAAA9TRkyBBNnTpV1dXVeu2113TzzTerVatWmjBhQsh5NTU1ys7OTsk127Rpk5L3qKysTEFrkAqMtCTDOJKMnJqqpm4JAABo4YzjyNlf0WRfxnHq1d6cnBx16tRJ3bt314033qhBgwbp1VdfDUzpeuihh9SlSxf16tVLkrRp0yaNGDFCpaWlatOmjYYPH64NGzYE3s+2bd1xxx0qLS1V27Zt9fOf/1zGmJBrhk8Pq66u1t13362ysjLl5OSoZ8+eeuaZZ7RhwwYNHDhQktS6dWtZlqWxY8dKipwetmfPHo0ePVqtW7dWfn6+hg4dqrVr1waenzZtmkpLSzV//nz17t1bhYWFGjJkiLZs2VKvzwvRMdKSQN03gZEsj1RV3dTNAQAALZw5UKndv/ppk12/zW+elFVY1ODX5+XladeuXZKkhQsXqri4WAsWLJAk1dbWavDgwRowYICWLFmirKwsPfjggxoyZIg+/vhjZWdn6/HHH9e0adP07LPPqnfv3nr88cf1yiuv6Pzzz495zdGjR2v58uX6wx/+oFNPPVXr16/Xzp07VVZWppdfflnf//73tWbNGhUXFysvLy/qe4wdO1Zr167Vq6++quLiYt199926+OKL9dlnn6lVq1aSpAMHDuixxx7T3/72N3k8Hl1zzTW68847NWPGjAZ/XqhDaEmGkeT1ytQQWgAAABrCGKOFCxdq/vz5uuWWW7Rjxw4VFBTo6aefDkwLmz59uhzH0dNPPy3LsiRJU6dOVWlpqRYvXqyLLrpITz31lCZMmKArrrhCkjRlyhTNnz8/5nW//PJLzZ49WwsWLNCgQYMkScccc0zgef9Usg4dOqi0tDTqe/jDyrJly3TWWWdJkmbMmKGysjLNmTNHP/jBDyTVha4pU6bo2GOPlSSNHz9eDzzwQEM/MgQhtCRyaLjRysqSqWZ6GAAAQH3MmzdPhYWFqq2tleM4uvrqq3X//ffr5ptv1imnnBJSx/LRRx/pq6++UlFR6EhOVVWV1q1bp3379mnLli3q379/4LmsrCz17ds3YoqY36pVq+T1enXuuec2+B4+//xzZWVlhVy3bdu26tWrlz7//PPAsfz8/EBgkaTOnTtr+/btDb4uDiO0JGLMoZGWLBmfT8b2yfLysQEAACRj4MCBmjx5srKzs9WlSxdlZR3+OaqgoCDk3P379+uMM86IOp2qffv2Dbp+rOlebvBPE/OzLCtmmEL98NN3QodqWrxeGduWfLZEaAEAAE3Eyi9Qm9882aTXr4+CggL17NkzqXNPP/10zZo1Sx06dFBxcXHUczp37qwVK1bonHPOkST5fD69//77Ov3006Oef8opp8hxHL355puB6WHB/CM9tm3HbFfv3r3l8/m0YsWKwPSwXbt2ac2aNTrxxBOTujc0DquHJWLq/mF5syTHZq8WAADQpCyPR57Coib7sjzu/fg4atQotWvXTsOHD9eSJUu0fv16LV68WLfeeqs2b94sSbrtttv0yCOPaM6cOfriiy900003ReyxEqxHjx4aM2aMrr/+es2ZMyfwnrNnz5Ykde/eXZZlad68edqxY4f2798f8R7HHXechg8frnHjxmnp0qX66KOPdM011+ioo47S8OHDXfksEIrQkgwjKcsr2bZEaAEAAHBFfn6+3nrrLXXr1k1XXHGFevfurRtuuEFVVVWBkZef/exnuvbaazVmzBgNGDBARUVFuvzyy+O+7+TJk3XllVfqpptu0gknnKBx48YF9mA56qijNHHiRN1zzz3q2LGjxo8fH/U9pk6dqjPOOEOXXHKJBgwYIGOMXnvttYgpYXCHZdI80a68vFwlJSXat29fzGG/5sTU1urgov+TPJbM/grlnn2+vG0bNqeyKfmHNPv37x8ylxRHNvo1M9GvmYu+zUxu92tVVZXWr1+vo48+Wrm5uSl/f0RnjFFlZaUKCgoCK5mhYeL9HU42GzDSktChTOfxyBgjU8tICwAAAJBOhJZEAptLWrIkyfY1cYMAAACAloXQkohRYKk6I8n4CC0AAABAOhFaEgkq+bEkielhAAAAQFoRWhIw/poWy5KxJKe2pmkbBAAAALQwhJZkHBptsTxeqaa6iRsDAAAAtCyElkSMqfuyLMmbJVNFaAEAAADSidCSiAn8Q/J6ZaqrmrI1AAAAQItDaEkoqBDf65Xx1cg4ThO2BwAAAGhZCC2JmMOF+PJ6JZ/DXi0AAKBZMNXVcir3p+XLVDffKfJjx47VZZddFnh83nnn6fbbb2/Uew4cOFB333134xqWwOLFi2VZlvbu3evqddxmWZbmzJnj6jWyXH33TGBM3d6SsmS8Xpnqaplan6xW2U3dMgAA0IKZ6mpVrVgiZ//+tFzPU1io3P5ny8rJSer8sWPH6rnnnpMktWrVSt26ddPo0aP1i1/8QllZ7v4I+ve//12tWrVK6tzFixdr4MCB2rNnj0pLSwPHX375ZdXUuLtq7FlnnaUtW7aopKQk6deMHTtWe/fudT0kNDeElmQYI1l1q4cZx5Z87NUCAACalvHVytm/X1Z2tqzs5IJEg69VUy1n/34ZX23SoUWShgwZoqlTp6q6ulqvvfaabr75ZrVq1UoTJkyIOLempkbZ2an5pXCbNm1S8h6VlZUpaE1s2dnZ6tSpk6vXiCWVn3c6MD0sKUbSoelhtiPD9DAAANBMWNk5snJz3f1qYCjKyclRp06d1L17d914440aNGiQXn31VUmHp3Q99NBD6tKli3r16iVJ2rRpk0aMGKHS0lK1adNGw4cP14YNGwLvadu27rjjDpWWlqpt27b6+c9/LhO0GbgUOT2surpad999t8rKypSTk6OePXvqmWee0YYNGzRw4EBJUuvWrWVZlsaOHSspcnrYnj17NHr0aLVu3Vr5+fkaOnSo1q5dG3h+2rRpKi0t1fz589W7d28VFhZqyJAh2rJlS8zPJ3x6WKL3uP/++/Xcc89p7ty5sixLlmVp8eLFSX1u0T7vX/ziF+rfv39Eu0499VQ98MADkqSVK1fqwgsvVLt27VRSUqJzzz1XH3zwQcx7cguhJZFATYskj7eunqWW0AIAAFBfeXl5IVOuFi5cqDVr1mjBggWaN2+eamtrNXjwYBUVFWnJkiVatmxZ4Ad3/+sef/xxTZs2Tc8++6yWLl2q3bt365VXXol73dGjR+v555/XH/7wB33++ef6y1/+osLCQpWVlenll1+WJK1Zs0ZbtmzRpEmTor7H2LFj9d577+nVV1/V8uXLZYzRxRdfrNrawzNwDhw4oMcee0x/+9vf9NZbb2njxo2688476/UZxXuPO++8UyNGjAgEmS1btuiss85K6nOL9nmPGjVK7777rtatWxc459NPP9XHH3+sq6++WpJUUVGhMWPGaOnSpXrnnXd03HHH6eKLL1ZFRUW97quxmB6WSCC516VZIyPD9DAAAICkGWO0cOFCzZ8/X7fcckvgeEFBgZ5++unANKXp06fLcRw9/fTTsixLkjR16lSVlpZq8eLFuuiii/TUU09pwoQJuuKKKyRJU6ZM0fz582Ne+8svv9Ts2bO1YMECDRo0SJJ0zDHHBJ73TyXr0KFDSE1LsLVr1+rVV1/VsmXLdNZZZ0mSZsyYobKyMs2ZM0c/+MEPJEm1tbWaMmWKjj32WEnS+PHjAyMWyYr3HoWFhcrLy1N1dXXItLJkPjcp8vOW6kZVZs6cqV/96leB++rfv7969uwpSTr//PND2vff//3fKi0t1ZtvvqlLLrmkXvfWGIy0JHKoEP8wi5oWAACAJMybN0+FhYXKzc3V0KFDddVVV+n+++8PPH/KKaeE/AD90Ucf6auvvlJRUZEKCwtVWFioNm3aqKqqSuvWrdO+ffu0ZcuWkClNWVlZ6tu3b8w2rFq1Sl6vV+eee26D7+Pzzz9XVlZWyHXbtm2rXr166fPPPw8cy8/PD4QNSercubO2b99er2s15D0SfW5+4Z+3JI0aNUozZ86UVBcun3/+eY0aNSrw/LZt2zRu3Dgdd9xxKikpUXFxsfbv36+NGzfW674ai5GWhIwOLR8mqe5fxrabskEAAABHhIEDB2ry5MnKzs5Wly5dIlYNKygoCHm8f/9+nXHGGZoxY0bEe7Vv375BbcjLy2vQ6xoifMUyy7Ii6m3ceI9kP7fwz1uSRo4cqbvvvlsffPCBDh48qE2bNumqq64KPD9mzBjt2rVLkyZNUvfu3ZWTk6MBAwa4vrJaOEJLMsyhQnxJRkaqZaQFAAAgkYKCgsA0o2ScfvrpmjVrljp06KDi4uKo53Tu3FkrVqzQOeecI0ny+Xx6//33dfrpp0c9/5RTTpHjOHrzzTcD08OC+Uce7Di/lO7du7d8Pp9WrFgRmB62a9curVmzRieeeGLS95cK2dnZEW1N5nOLpWvXrjr33HM1Y8YMHTx4UBdeeKE6dOgQeH7ZsmX685//rIsvvlhSXcH/zp07G38j9cT0sATMoYEW/0iLPF451VVN2CIAAIDMNGrUKLVr107Dhw/XkiVLtH79ei1evFi33nqrNm/eLEm67bbb9Mgjj2jOnDn64osvdNNNN8XdnLFHjx4aM2aMrr/+es2ZMyfwnrNnz5Ykde/eXZZlad68edqxY4f2R9n35rjjjtPw4cM1btw4LV26VB999JGuueYaHXXUURo+fLgrn0W8+/n444+1Zs0a7dy5U7W1tUl9bvGMGjVKL7zwgl588cWQqWFS3b3/7W9/0+eff64VK1Zo1KhRaR298iO0JBJUiC9JljdLhtACAACaCVNTLVNV5e5XTXVa7iU/P19vvfWWunXrpiuuuEK9e/fWDTfcoKqqqsAIws9+9jNde+21GjNmjAYMGKCioiJdfvnlcd938uTJuvLKK3XTTTfphBNO0Lhx4wJ7sBx11FGaOHGi7rnnHnXs2FHjx4+P+h5Tp07VGWecoUsuuUQDBgyQMUavvfZa0ptYpsq4cePUq1cv9e3bV+3bt9eyZcuS+tziufLKK7Vr1y4dOHBAl112WchzzzzzjPbs2aPTTz9d1157rW699daQkZh0sUx9J9o1Unl5uUpKSrRv3756D181BXvPLlUtWSRP+/ayvFmyd+2Up7hEed85r6mbVi/+Ic3+/fu7vgst0od+zUz0a+aibzOT2/1aVVWl9evX6+ijj1Zubm7guKmuVtWKJXKijAy4wVNYqNz+Z9drc8kjmTFGlZWVKigoCKzIhYaJ9XdYSj4b8F/MRALzww7xemTSXHgEAAAQzsrJUW7/s9O2FYOV1arFBBY0P4SWRIxCCvEtb5Zk18r4fLL4LRkAAGhCVk4OQQItAjUtCZm60OIfFfR6ZWyHDSYBAACANCG0JM0/0uKVfL66LwAAAACuI7Qk4l+nIGikRY4tQ2gBAABplOa1k4CUScXfXUJLIsYELXssyeOVbFtiehgAAEgD/5K6Bw4caOKWAA3j/7vbmOWhqSRPJJBXDk0P83hkjGGkBQAApIXX61Vpaam2b98uqW4vE5bgdZ8xRtXV1fJ6vXzeDWSM0YEDB7R9+3aVlpbK6/U2+L0ILQkZGZmQv6yWMVItIy0AACA9OnXqJEmB4AL3GWNUU1Oj7OxsQksjlZaWBv4ONxShJYFoc/CMZcnYjLQAAID0sCxLnTt3VocOHVTLL07Twufz6eOPP1bPnj3ZDLYRWrVq1agRFj96IJGgPVr8LInVwwAAQNp5vd6U/ACIxHw+n4wxys3NJbQ0AxTiJyV0tMV4LDnV1U3UFgAAAKBlIbQkJWykxeOVaqqaqC0AAABAy0JoSSTautJerwwjLQAAAEBaEFoSMUbh08MILQAAAED6EFoSilKI782Sqa2Vse2maRIAAADQghBakhCxMrfXKzm2xLLHAAAAgOsILYlEq2nxeGRsnwzLHgMAAACuI7QkYhQRXCxvlmQ7Eps7AQAAAK4jtCRijGSFTRDzeiVGWgAAAIC0ILQkZMLXDqubHkZNCwAAAJAWhJZEotS0WJYlyZLxMT0MAAAAcBuhJYFodfiB52oZaQEAAADcRmhpKGPqlj0GAAAA4CpCS0LRh1qs8OJ8AAAAAK4gtCRiTOw5Yo6T3rYAAAAALRChJZFoSx4HPwcAAADAVYSWBoqyEDIAAAAAFxBaEjFGMatXGGkBAAAAXEdoSSRmMLEILQAAAEAaEFqSES2cWJKhEB8AAABwHaElkZiF+JaMIbQAAAAAbiO0NJRlSQ7TwwAAAAC3NSq0PPLII7IsS7fffnuKmtMMGRN7pTDHTm9bAAAAgBaowaFl5cqV+stf/qI+ffqksj3NT6xie4tCfAAAACAdGhRa9u/fr1GjRul//ud/1Lp161S3qRmKrGmxJImaFgAAAMB1DQotN998s4YNG6ZBgwaluj3NjpGJllkky0NNCwAAAJAGWfV9wQsvvKAPPvhAK1euTOr86upqVVdXBx6Xl5dLknw+n3w+X30vn3a2z5YtyQqbCmZLsnz2EXEPkmTbtowxsm3qcDIJ/ZqZ6NfMRd9mJvo1M9Gv6ZHsz9L1Ci2bNm3SbbfdpgULFig3Nzep1zz88MOaOHFixPH33ntPBQUF9bl8k3Aq98tRtqzd5WFPeKWdu+VdsaJpGlZPjuOooqJC7777rjweFo3LFPRrZqJfMxd9m5no18xEv6ZHZWVlUudZxiRfTT5nzhxdfvnl8nq9gWO2bcuyLHk8HlVXV4c8J0UfaSkrK9OuXbtUXFyc7KWbTM3nn6j232uU1emokOP2ju3ytu+onDP6N1HL6se2ba1cuVL9+vWL6CMcuejXzES/Zi76NjPRr5mJfk2P8vJytW3bVvv27YubDeo10nLBBRdo9erVIceuu+46nXDCCbr77rujdmhOTo5ycnIiL5yVpayses9OSzvbIzmy5A3fYNJjySsdEffgZ1mWvF7vEdVmJEa/Zib6NXPRt5mJfs1M9Kv7kv1s69UDRUVFOvnkk0OOFRQUqG3bthHHM0bMYnuL1cMAAACANGCCXjLCR1mkQ5mF0AIAAAC4rdFjXYsXL05BM5oxY6KueFyXWljyGAAAAHAbIy2JxJoCZjE9DAAAAEgHQksijok5PUxMDwMAAABcR2hJxBgp6gQxi9ACAAAApAGhJRFjYmYWSloAAAAA9xFaEjDGUbTUYlkealoAAACANCC0NIqRYbgFAAAAcBWhJZGY08MsyRFzxAAAAACXEVoScWIV4kuSYYoYAAAA4DJCSyLGiT3SIiMx0AIAAAC4itCSSKwlj5keBgAAAKQFoSWRWDUtdU+KoRYAAADAXYSWBOoiSfSRFmMMIy0AAACAywgtiSQKJQ6hBQAAAHAToSWRuIX4kmF6GAAAAOAqQksisZY8tizJcZgeBgAAALiM0JJIokJ8MgsAAADgKkJLIsZRzJGWuhPS2RoAAACgxSG0JCNqZrHq8orjpL05AAAAQEtCaEmgbnZY1Ep86lkAAACANCC0JGJijKRYUl1NC8EFAAAAcBOhJQ5jDu14b8UaaRHTwwAAAACXEVriMXFWB7Ms1T3JSAsAAADgJkJLIkbRR1osSY5hdhgAAADgMkJLPCbeSIoVdA4AAAAAtxBa4ooTSPyF+EwPAwAAAFxFaInHHPpHzEJ8IzmEFgAAAMBNhJZ4EhXix50+BgAAACAVCC1xmfiF+BI1LQAAAIDLCC3xJCrEZ3oYAAAA4DpCSwLGRK9psSxLhkJ8AAAAwHWElngS5hGrLtQAAAAAcA2hJZ5kAgmhBQAAAHAVoSUO4x9qibrkMQAAAIB0ILTEY0yCkRQjOU7amgMAAAC0RISWRGIU4texRCE+AAAA4C5CSzwm8I+oLImaFgAAAMBlhJa4kqhpIbMAAAAAriK0xHNoFMVS9NBSNxBDTQsAAADgJkJLPMbUJZN4i4cxPQwAAABwFaElkQShhMwCAAAAuIvQEo8xSmKoJU2NAQAAAFomQksy4mUW9mkBAAAAXEVoicdf0xI3tQAAAABwE6ElHn/BSozMwj4tAAAAgPsILXH5a1rinUJoAQAAANxEaImjrg4/QSE+oQUAAABwFaElngT7tBgZQgsAAADgMkJLXP5AwkgLAAAA0FQILQmZ2JnFsupGWwAAAAC4htASjzEJalosGdtOZ4sAAACAFofQEk+iQRRLksNICwAAAOAmQktch0ZaYpa0WJKcNLYHAAAAaHkILfGYBIX4lsVICwAAAOAyQksyYmUWSTKMtAAAAABuIrTEk6gQ37Ikm9ACAAAAuInQEs+hmV+WFXv1MPZpAQAAANxFaInLxN+HxSK0AAAAAG4jtMRhEgUSSzLUtAAAAACuIrTEE3djSdU9R00LAAAA4CpCS2NYYvUwAAAAwGWElsagpgUAAABwHaElHmOkeIX4shLXvQAAAABoFEJLPEkU4sux09IUAAAAoKUitMSVRCG+w0gLAAAA4CZCSzzGxI0sdYX4hBYAAADATYSWeBLlEcti9TAAAADAZYSWeBKMoliyZJgeBgAAALiK0JJA3EhiWZIcVhADAAAAXERoiSvRksf+UwgtAAAAgFsILfGYBKuHWRahBQAAAHAZoSWOhFnEspTUaAwAAACABiO0xJVEGDGGzAIAAAC4iNDSGIHpYSx7DAAAALilXqFl8uTJ6tOnj4qLi1VcXKwBAwbo9ddfd6ttTc+Y+HPE/NPDqGkBAAAAXFOv0NK1a1c98sgjev/99/Xee+/p/PPP1/Dhw/Xpp5+61b6mZcyhYBLvFKaHAQAAAG7Kqs/Jl156acjjhx56SJMnT9Y777yjk046KaUNOyL4Aw0jLQAAAIBr6hVagtm2rRdffFGVlZUaMGBAzPOqq6tVXV0deFxeXi5J8vl88vl8Db18Wti2LUeSHSOUOMbIOI5qbZ88R8C9GGNk23ZTNwUpRL9mJvo1c9G3mYl+zUz0a3okmwcsU8/t3FevXq0BAwaoqqpKhYWFmjlzpi6++OKY599///2aOHFixPH58+eroKCgPpdOO+fgAZn9FVJObowTbMl25CltLcvrTW/j6slxHFVUVKioqEgeD+svZAr6NTPRr5mLvs1M9Gtmol/To7KyUoMHD9a+fftUXFwc87x6h5aamhpt3LhR+/bt00svvaSnn35ab775pk488cSo50cbaSkrK9OuXbviNqw58P17rWo+/UjeLl2jPu9UHZSprFTu2efLk5ef5tbVj23bWrlypfr16ydvMw9YSB79mpno18xF32Ym+jUz0a/pUV5errZt2yYMLfWeHpadna2ePXtKks444wytXLlSkyZN0l/+8peo5+fk5CgnJyfywllZyspq8Oy0tDAej7ySvDGK8T2WR46kLI9HnmZ+L5JkWZa8Xm+z/9xRP/RrZqJfMxd9m5no18xEv7ov2c+20WNdjuOEjKRkGhNvaTBLiZdFBgAAANAo9YqNEyZM0NChQ9WtWzdVVFRo5syZWrx4sebPn+9W+5pWojDCPi0AAACA6+oVWrZv367Ro0dry5YtKikpUZ8+fTR//nxdeOGFbrWvaSUMI1bdHi2EFgAAAMA19QotzzzzjFvtaJbqpobF2VzSqjuL0AIAAAC4h/Xb4nFM3MxSN9JipHh1LwAAAAAahdAShzGOEo60GMk4hBYAAADALYSWROKNtASWQia0AAAAAG4htMSTcATFkoxDTQsAAADgIkJLXMlMDzMMtAAAAAAuIrTEk2ikxT89zDjutwUAAABooQgt8RgTVLcSjXX4PAAAAACuILTEY5wEdfhW3QpjZBYAAADANYSWeJIOI6QWAAAAwC2ElngSTg+TJEvGoaYFAAAAcAuhJR5jFH+jFgAAAABuI7TEY0wSmcVQiA8AAAC4iNAShzEJ9mmR6p4ntAAAAACuIbTEk8RIS93ThBYAAADALYSWhJKoaUm0CSUAAACABiO0xOMkU9MipocBAAAALiK0xGMSL2Vsgv4JAAAAIPUILfEkueSxYXoYAAAA4BpCSxwmqSWPJUZaAAAAAPcQWuJJZqTFsE8LAAAA4CZCSzzGSWKbFvZpAQAAANxEaEkofmpJavYYAAAAgAbLauoGNGfm4EFVL39Lzv5yterZS9mn95dlRYkpTuJVxgAAAAA0DCMtcdR89rF8mzbIVJSr5sOVsr/eFP1EpocBAAAAriG0xOHs2hHy2N6yOeIcw8phAAAAgKsILXEYO7QQ39m3N/qJTA8DAAAAXENoicfxKTi1OOX7opxkMdoCAAAAuIjQEk/YCIpTvq9uw8lgliSH0AIAAAC4hdASgzEmctqXr1bm4IHQY5YlY5geBgAAALiF0BKLMXU1LWEi61osaloAAAAAFxFaYok20iLJRKtrIbQAAAAAriG0xGSkKNO+IorxLYt9WgAAAAAXEVriiTKC4pTvDXlsSVHDDQAAAIDUILTEYiTj2BGHI0daPKweBgAAALiI0BJLjJoWp3xv6LLHlqhpAQAAAFxEaIkpemhRba1M1cGgAyx5DAAAALiJ0BJLjJEWSTLByx5TiA8AAAC4itASh3GculASJqQY3xI1LQAAAICLCC0xGNuJOYISWoxvSVEK9gEAAACkBqElFp8v5lMhocVSaGE+AAAAgJQitMRgnDihJbimRdS0AAAAAG4itMRgEoy0BEZXLIvNJQEAAAAXEVpi8cWpU6mtObzsMfu0AAAAAK4itMRixx5pkSQTqGuxCC0AAACAiwgtMZgEK4L5i/Ety6KkBQAAAHARoSWWQzUtkbu01Ans1UJNCwAAAOAqQksM8QrxpfAVxAzLHgMAAAAuIbTEYic3PUyWJTli2WMAAADAJYSWWGyfJKNYE8Sc8r1BoyuGKWIAAACASwgtMQSmh8UqaqmpkamuqhtpkanLNwAAAABSjtASi88XEkSsnBzJCv24zL59TA8DAAAAXEZoiSV875WsVvIUFYee4l9BTEYMtQAAAADuILTEYOza0AOWR1ZJScghp7xupMUYw0gLAAAA4BJCSwzGF7p6mOXxyFNcGnLs8EiLJIfQAgAAALiB0BJL+JLHUUPLoZoWSYbpYQAAAIArCC0xGCdsc0mPR57i0OlhpnzvoUJ8h+lhAAAAgEsILbHYYYX4Ho88JaUhh0x1tUx1teqWPCa0AAAAAG4gtMRih460WF6vrMKiyGWP91cc+gOhBQAAAHADoSUGYzsKWcbY46krxi8qCjnPqdh3aMVjQgsAAADgBkJLLD7/ksd1hfby1H1UVnhdS2WlZMKmkgEAAABIGUJLLFFWD5Mkq1VOyOHAfi6MtAAAAACuILTEYPyh5dBAi+Xx1v0hyxt6os+umx7mMNoCAAAAuIHQEovtCylpCYy0eLNCz3NsHSpqSVfLAAAAgBaF0BJLlCWPJUlZoaHF2LbkGGaHAQAAAC4htMRgnNCaFv/0MCsstMh3aGlkUgsAAADgCkJLDMYXuk9LYKTFG1rTYhz/0siEFgAAAMANhJZY7OihJaKmxfbVjbI4hBYAAADADYSWWMJXA/OPsISFFuOzD00NI7QAAAAAbiC0xGDs8JqWQyMt4TUtNjUtAAAAgJsILbHEqmmJFlqYHgYAAAC4htASS9jqYbEK8WXbMhTiAwAAAK4htMQQPj0sViF+3XmWDNPDAAAAAFfUK7Q8/PDD6tevn4qKitShQwdddtllWrNmjVtta1rh+7RYMaaH+Wrr/k1oAQAAAFxRr9Dy5ptv6uabb9Y777yjBQsWqLa2VhdddJEqKyvdal/TscNWD/NvLhlt9TAAAAAArslKfMphb7zxRsjjadOmqUOHDnr//fd1zjnnpLRhTS1iepjXP9ISXtPikzFO5BLJAAAAAFKiXqEl3L59+yRJbdq0iXlOdXW1qqurA4/Ly8slST6fT77wFbqaEePzhZTWO5ZHtjFyPN6IknvHkWzbjlxxrBmxbVvGmLp2ImPQr5mJfs1c9G1mol8zE/2aHsnmAcs0sILccRx973vf0969e7V06dKY591///2aOHFixPH58+eroKCgIZdOi/ZzZip751bpUC3Lrn7fVWWP4+Q9UKmj/jE75NxNQ74vtWsvT05uUzQ1KY7jqKKiQkVFRfJ4WH8hU9CvmYl+zVz0bWaiXzMT/ZoelZWVGjx4sPbt26fi4uKY5zV4pOXmm2/WJ598EjewSNKECRN0xx13BB6Xl5errKxMffv2jduwprbv/16RzxhZreqmg7UuKVRWm2KZ/GwdCJsi1sc+qNxjj5H3qLKmaGpSbNvWypUr1a9fP3nDl23GEYt+zUz0a+aibzMT/ZqZ6Nf08M/CSqRBoWX8+PGaN2+e3nrrLXXt2jXuuTk5OcrJyYm8cFaWssJX4mpGLMeRJck69Njj8cprWTKtsgLHAufatrweq1nfjyRZliWv19vs24n6oV8zE/2auejbzES/Zib61X3Jfrb16gFjjG655Ra98sorWrx4sY4++ugGNe5IYJLdXFKSbB9LHgMAAAAuqVdoufnmmzVz5kzNnTtXRUVF2rp1qySppKREeXl5rjSwyYQVXVmHwopleeqWPw4ONbZDZgEAAABcUq+qosmTJ2vfvn0677zz1Llz58DXrFmz3Gpf0wlfKSK4ACuspqVuVIbUAgAAALih3tPDWoqY08NUt8GkUc3h52ybfVoAAAAAl7B+Wyx2aAixrKCPyhuW9QgsAAAAgGsILbGEj7QEFeBb4ascUIgPAAAAuIbQEoOJW9MSHlpsQgsAAADgEkJLFMaYyClfITUtYYX4hBYAAADANYSWaMJHWSRZnqCgElbTYhxCCwAAAOAWQks0Pl/kseCRFqaHAQAAAGlDaIkiYrljKbSmJcrqYcawghgAAADgBkJLNFGmh8XdXNK2ZVj2GAAAAHAFoSWKupXDQqd7Bde0WBEjLbbkMD0MAAAAcAOhJRr/SIsVdCzukseOJEZaAAAAADcQWqKx7fCBlrAlj6OsHsZICwAAAOAKQksUCQvxw2pa6kIOIy0AAACAGwgt0QSWPD40P8zyyLIOzxWLrGlxDk0RAwAAAJBqhJYoIgrxPWEfU3hoYZ8WAAAAwDWElmjCp4d5Qz8mK8qSx4QWAAAAwB2Elmh8vpCBluDljiVF2VzSjl4HAwAAAKDRCC1RJJweFrHkMauHAQAAAG4htERjh42ahIWWqEses3oYAAAA4ApCSxTG9oU8tpIZaaGmBQAAAHAFoSWaiJGW0JoWyxtZiG8ILQAAAIArCC3R2HZISUvCJY8dO2J0BgAAAEBqEFqiqAsgsQvxrfDpYUZBG1ICAAAASCVCSzT+6WHWoX8lqmmRCC0AAACASwgtUZhATcuh1BJR0xIWWizJ1Na63zAAAACgBSK0RJOwpiVss0lZMr4at1sFAAAAtEiElmjCd7ePmB4WHlok+XysIAYAAAC4gNASRd30sMMBJLymxbI8EVPGDHu1AAAAAK4gtESTYJ8WSZGjLT4foQUAAABwAaEliog9V8KnhymyGD9imWQAAAAAKUFoiSZipCXKxxS+gpjPR2YBAAAAXEBoiSZs9TDLG2WkJXyvFp8tGcflhgEAAAAtD6ElCuMLm+oVraYleNljy6qbHkZNCwAAAJByhJZoEi15rMiRFsP0MAAAAMAVhJYoTFhNS/iSx5Iia1pY8hgAAABwBaElGp+vbtDEOvQ4yvSwyJoWnww1LQAAAEDKEVqicfyjJodSS9SRlvDNJZkeBgAAALiB0BJF+PSwpJY8Zp8WAAAAwBWElijCN5eMVtMSPD3M0qGgQ00LAAAAkHKElmh8oaEl6pLHUWpaCC0AAABA6hFaovElseRx8PQwy2L1MAAAAMAlhJYojBM+0pJMIT6hBQAAAHADoSWaJPZpscJCi2ymhwEAAABuILREk8zqYVmtoryG0AIAAACkGqElioglj8NHVSRZWZHTw4xDaAEAAABSjdASTcTqYcns08JICwAAAOAGQksU4SMt0WpaIkZfqGkBAAAAXEFoica/Ephl1T2OurlkaE0Lq4cBAAAA7iC0RGMns7lk0DHLYqQFAAAAcAmhJYqIQvxEm0tKku0QWgAAAAAXEFqiSaqmJTS0GMfH6mEAAACACwgtURgnfKQl2pLHYSMtRpJd616jAAAAgBaK0BLGGJPc5pJR9m5xampcahUAAADQchFawoWPsijG9LDwkRZLUi0jLQAAAECqEVrChY+ySFFHVSIK8Y1kfIQWAAAAINUILWGM73Bosfx/iDrSEmUZZKaHAQAAAClHaAnnRNkkMtqSx5YnokDf1BJaAAAAgFQjtIQL2Viybqwlak2LFDltjJoWAAAAIOUILWEOTw8zh+eHxQgt4cseG0ILAAAAkHKElnDRCvFjjrSEFeNTiA8AAACkHKElTGBjyeCyliRHWpwaQgsAAACQaoSWcCE1LZIsq67oPprgmhZLjLQAAAAALiC0hAufHhZrapgkKzi0GMn4WD0MAAAASDVCSxhTj9CirFahjynEBwAAAFKO0BLODt2nxfJE2UTS/1xW+D4thBYAAAAg1QgtYUx4TUu8kZaI1cN80c8DAAAA0GCElnD1mR7mZZ8WAAAAwG2ElnARoSXe9DD2aQEAAADcRmgJc7gQv66uxYo70kJNCwAAAOA2Qku4kJEWK/6Sx2EjLYaRFgAAACDlCC1hQgrxLSVf02KJJY8BAAAAFxBawoUveeyNU9MSPj0sfOUxAAAAAI1GaAnnnx7mzy1Jby5pydTWuNUqAAAAoMWqd2h56623dOmll6pLly6yLEtz5sxxoVlNx9RryeOwUZhaRloAAACAVKt3aKmsrNSpp56qP/3pT260p+lFbC6Z/JLHFOIDAAAAqZeV+JRQQ4cO1dChQ91oS7MQPtJSnyWP5WOkBQAAAEg1alrC1WN6mBVS08JICwAAAOCGeo+01Fd1dbWqq6sDj8vLyyVJPp9PvmY4MmHX1sgYE6jDN5ZHdtBqYsEcr1fBz5ja5nlPkmTbtowxssNDGY5o9Gtmol8zF32bmejXzES/pkeyPzu7HloefvhhTZw4MeL4e++9p4KCArcvX28l//63CquqJcuSfLYqa2q1a3d51HOz9x9UJ5//L7IlU1WlL1esSF9j68FxHFVUVOjdd9+VJ96UNxxR6NfMRL9mLvo2M9GvmYl+TY/KysqkznM9tEyYMEF33HFH4HF5ebnKysrUt29fFRcXu335eju49T+q+vITOZW1slp5VZyfpx5torfTMdU6mHWorqXWJ1nSmWeeKcuy0tji5Ni2rZUrV6pfv37yxtl7BkcW+jUz0a+Zi77NTPRrZqJf08M/CysR10NLTk6OcnJyIi+claWsLNcvX2+WjCwp8OXxeuWNEUKsVtnyP2OsutdmKXJVsebCsix5vd5m+bmj4ejXzES/Zi76NjPRr5mJfnVfsp9tvXtg//79+uqrrwKP169fr1WrVqlNmzbq1q1bfd+u+fE1Yp8WY2Rqa2W1ahX9fAAAAAD1Vu/Q8t5772ngwIGBx/6pX2PGjNG0adNS1rAm49QjtERLhqwgBgAAAKRUvUPLeeedJxNjNa1MUJ99Wqwo8xtNLaEFAAAASCWWQggXvuxavMKrKCMtprYmxQ0CAAAAWjZCS5jwkZa4m0tantDnjZEILQAAAEBKEVrCObYUvGVkonW5vaGjLUwPAwAAAFKL0BLOthVcsmN54q/LHVje2LJkJBkK8QEAAICUIrSEMeE1LfUZaTFGqmF6GAAAAJBKhJZw9VnyWJEriDmEFgAAACClCC1hIpc8jj89LGIFMQrxAQAAgJQitISr5/Sw8JEWU1Od6hYBAAAALRqhJVw9ljyWJGW1CnnI6mEAAABAahFawpj61rRkhU0fo6YFAAAASClCSzg7dJ8WK9FIiydsnxaWPAYAAABSitASJrwQXwkL8Q89b1kseQwAAAC4gNASzraDB1qSmB4WXtNCaAEAAABSidASzq7v5pLhq4cRWgAAAIBUIrSEqe8+LVYWNS0AAACAmwgt4eq75HHYSAubSwIAAACpRWgJYoyJLMT3Jtpcsm6kxfK/B/u0AAAAAClFaAnmOKqrwjeHU0jCzSXDpodR0wIAAACkFKElWPjGkkpc0yKvP7RYkozESAsAAACQUoSWIIGpYfVa8vhQaDk0MsP0MAAAACC1CC3BfL7IY/UpxDeSWD0MAAAASClCS7BAEb6Rf+ik3kses3oYAAAAkFKEliAmSk1L4pGW8NDCSAsAAACQSoSWYNGmhyW55LEfoQUAAABILUJLkIYU4isrbPqYr1bGmOjnAgAAAKg3QkswJ7imRZIsWVai1cNahR4wJqg2BgAAAEBjEVqChU8PSzA1rO4c/0jLoX1aRDE+AAAAkEqEliAmfIQk0dQw6XBosYKOUdcCAAAApAyhJVhYTYuVRGgJWfLYSDKGYnwAAAAghQgtQSKWPE6wR4skKTi0HBptYXoYAAAAkDqElmB2WCF+MiMtlufwef76fUILAAAAkDKElmDBNS2WpAQrhwV4Q0dkTLT9XgAAAAA0CKElSGghviUrmdXDFLbssTFMDwMAAABSiNASzD40QuKf5pVMTYsUMtJijJFqCC0AAABAqhBagjRoyWMFjbRYkhxHpro6tQ0DAAAAWjBCS7AGFOJLkpVf4P+T5Niy9+xMedMAAACAlorQEswOLaBPZp8WSfIUFR16gSU5jpxdu1LdMgAAAKDFIrQEiZwellxNi1VUfOgPlowxsnduS3HLAAAAgJaL0BKsgTUtnsK60GJZHsk4cnbtSHXLAAAAgBaL0BLEBFYPq6tpSXbJY09gpKXutfbuXTKO40ILAQAAgJaH0BIsBdPDZCTV1sop35fatgEAAAAtFKElWHhosZJcPSwvL2KvFmc3K4gBAAAAqUBoCRJSiG9JSnJ6mGV55Cn0ryCmuhXECC0AAABAShBagh2qaTlU0ZL0ksdS8BQxT91eLbtZ9hgAAABIBUJLkMBIi/FvLplcTYt0eAUxWZaM48hmBTEAAAAgJQgtwRq45LEUvILYoQ0md25PYcMAAACAlovQEsy2/XPD6jRgephlWZKMbEILAAAAkBKEliAmbKSlPjUtwdPDZIycPbsj3g8AAABA/RFagjm2QoZa6jXSUnT4gZFk23L27UlZ0wAAAICWitASxPh8oQfqUYhv5eZJWVmH38s4clhBDAAAAGg0QkuwRhTiW5YVOkXMcWSzVwsAAADQaISWYE7Da1qk4L1aLMlmpAUAAABIBUJLEOMLH2lJfnqYFLrssXFsVhADAAAAUoDQEswOr2mp38cTvoIYoQUAAABoPEJLkJAlii0rEFqMbcsYJ+Hr/SuIWf5lj3fucKWdAAAAQEuSlfiUFsQ5tLmkMZJlyfLWTQ+zt2yWjORp116evPyYLw+eHiZj5JTvlfH5ZGXxMQMAAAANxUhLMF/kPi3GsWV5vMrq1kNm7x7ZO3fIGBP15YHQ4mfbcvbudq+9AAAAQAtAaAkSsYO9xyNTVSXl5in7pFOV0/fbsnJy5Hy9WaamJuL1Vk6ulJ19+P1Y9hgAAABoNOYtBYuy5LGpOihPfoGsvHxl5RfIU9JatV98qtqN6+Vp3UaegsKQ13gKi+Xs3nloipgjZxehBQAAAGgMRlqCGF/46mFemaoqedq2qyuul+QpLFL2aX3VqteJcvZETv0KqWuxbfZqAQAAABqJ0BIsyvQwOY68RSUhhy1vlrI6HSUrq5VMbeg0MStkrxZH9i5WEAMAAAAag9ASLGx6mGRkeb2yCooiTvWUtpanpFROeXnocfZqAQAAAFKK0BIkuBDfkmRqa6XcvIi6FUmyvF55j+oqc/BAyHFPxF4thBYAAACgMQgtwWy7bo+WQ0xtrTyFRbJyc6Oe7m3bXlZ2tkx1deBY8PQwycgp31cXfgAAAAA0CKElWHhNi69WnjZtY57uKWktT2lrORWHp4iF7NVyaKNKZw/F+AAAAEBDEVoOMY4jYxz/I0mWLFnyFhbHfI3l8SirS5lM1cHDx1pl1+3XEvS+7NUCAAAANByhxS9slMXISK2yZRVG1rME87ZtJysnp24TykNCpog57NUCAAAANAahxS84tBhJjiMrLz9qEX4wq7hUntZt5VTsCxwLWUHMsRlpAQAAABqB0HKICV/u2DjyFBfLys6J+zrLspTVpatMVZXMoSJ+/wpi7NUCAAAANB6hxc/nC31sjDzt2if1Um+bdrLy8gO1LSHTw4yRs4NljwEAAICGIrQcYsJXDpMlb1FJUq+1iorlbdtW5tAqYv4VxPx7tTDSAgAAADQcocUvfHqYZYUuXxyHZVnyduoq1VTLGCOrMGykpXxfyF4uAAAAAJLXoNDypz/9ST169FBubq769++vd999N9XtSjvjnx7m31zS45FVXJr0671t20l5+TIHKg/XtBx6H1N1QJVvzJVTW5O6BgMAAAAtRL1Dy6xZs3THHXfovvvu0wcffKBTTz1VgwcP1vbtR3jdRtj0MMvjkZWdnfTLPYVFyurURWbPbllZrWTl5de9T1aWZHl04JUXtPvnN2v/nFmyd2wLFO0DAAAAiC+rvi944oknNG7cOF133XWSpClTpugf//iHnn32Wd1zzz0pb6AbKiv3aus3a0OOebbvULZTKxlbjuXIauXVuqpt9Xpfc1SJqnYYWfu/VmHbYmVtrqh7wqu6EZe9O1Q170Xt++erUm6eTOvWclq3lmndRqa4WMrNlcnJkXKyZbJzpFatJCv4ClboBa2wxzFPtWQ7tvbu/Err13rk9da729FM2bZPe3euo18zDP2auejbzES/ZqaW0K+l7bqoTZujmroZSbFMPX7lX1NTo/z8fL300ku67LLLAsfHjBmjvXv3au7cuQnfo7y8XCUlJdq3b5+Ki5OrGUm1Lz5epIo/PxH1OSMjOY5qs7P0x5vPq/d7m6oqOQcPqMvual2+4EvlHawNO8FIsuoChyXJsmSFhxEX1G09Y8vj8abhakgX+jUz0a+Zi77NTPRrZmoJ/Vp78cU663s3Nmkbks0G9YqNO3fulG3b6tixY8jxjh076osvvoj6murqalUHFaGXl9etsOXz+eQLX2Y4TWzbUSCpGSf0yUMZzvF4GjaFKztbqqnRN+08evqqU9X7q106ec0Oddqxv+55/+iIceq+G2TJWHX/VtC/XGOC7h2Zg37NTPRr5qJvMxP9mpkytV8tj4xjmuzncb9kr+/6WNfDDz+siRMnRhx/7733VFBQ4Pblo9q1ZZ06+FcLi/G3cEtpriorKxv0/kaWZKQDHktvH99Obx/fTh13HdC3vtyhk9ftUl61rZB0YoIa4vJ3hUMtTUaiXzMT/Zq56NvMRL9mpoztV8vWtu3btGLFiiZtRrI/b9crtLRr105er1fbtoXWemzbtk2dOnWK+poJEybojjvuCDwuLy9XWVmZ+vbt22TTw9Z+WqO9b8Yust/VtlDLLjxJrQsbGqqMzP4KmdraQDF/TcdsrWhXpHf7lamwVio+aKtkf42KK6pVvL9aeQdrlF1rq1WtrVY1trJrbXntw6NAVsT3i4nzXLQmGdmOI6/Hk7FDnC2SkWzHltfjdX+UDulDv2Yu+jYz0a+ZqQX0a8fOndW/f/8mbYN/FlYi9Qot2dnZOuOMM7Rw4cJATYvjOFq4cKHGjx8f9TU5OTnKycmJvHBWlrKymqaoqfepA6X/Ghj3nEsaeQ1753ZVvbNUVmGhPPkFsvfskqqqlHXscco+/kRZ2ZGfiZt8Pp9WrFihM/v3b7LPHann79d+9GtGoV8zF32bmejXzES/pkeyn229e+COO+7QmDFj1LdvX5155pl66qmnVFlZGVhNDHU8bdsrq9vRqv3qCzn79srKzVPOaf3k7dpdloc9PQEAAIBk1Tu0XHXVVdqxY4d+/etfa+vWrfrWt76lN954I6I4v6WzLEutjukpe/tWWTk5yj6pj7xt2jV1swAAAIAjToPGusaPHx9zOhgO8xQWKeeM/rLy8uTJzWvq5gAAAABHJCbouczbuk1TNwEAAAA4olFcAQAAAKBZI7QAAAAAaNYILQAAAACaNUILAAAAgGaN0AIAAACgWSO0AAAAAGjWCC0AAAAAmjVCCwAAAIBmjdACAAAAoFkjtAAAAABo1ggtAAAAAJo1QgsAAACAZo3QAgAAAKBZI7QAAAAAaNYILQAAAACaNUILAAAAgGaN0AIAAACgWctK9wWNMZKk8vLydF+6RfP5fKqsrFR5ebmystLe7XAJ/ZqZ6NfMRd9mJvo1M9Gv6eHPBP6MEEvae6CiokKSVFZWlu5LAwAAAGiGKioqVFJSEvN5yySKNSnmOI6++eYbFRUVybKsdF66RSsvL1dZWZk2bdqk4uLipm4OUoR+zUz0a+aibzMT/ZqZ6Nf0MMaooqJCXbp0kccTu3Il7SMtHo9HXbt2TfdlcUhxcTHfeBmIfs1M9Gvmom8zE/2amehX98UbYfGjEB8AAABAs0ZoAQAAANCsEVpaiJycHN13333Kyclp6qYghejXzES/Zi76NjPRr5mJfm1e0l6IDwAAAAD1wUgLAAAAgGaN0AIAAACgWSO0AAAAAGjWCC0AAAAAmjVCyxHkrbfe0qWXXqouXbrIsizNmTMn5Plt27Zp7Nix6tKli/Lz8zVkyBCtXbs25JzzzjtPlmWFfP3kJz8JOWfjxo0aNmyY8vPz1aFDB911113y+Xxu316LlY5+/eijjzRy5EiVlZUpLy9PvXv31qRJk9Jxey1Wur5f/Xbt2qWuXbvKsizt3bvXpbtCOvt12rRp6tOnj3Jzc9WhQwfdfPPNbt5ai5euvl25cqUuuOAClZaWqnXr1ho8eLA++ugjt2+vxUpFv0rS8uXLdf7556ugoEDFxcU655xzdPDgwcDzu3fv1qhRo1RcXKzS0lLdcMMN2r9/v9u316IQWo4glZWVOvXUU/WnP/0p4jljjC677DL9+9//1ty5c/Xhhx+qe/fuGjRokCorK0POHTdunLZs2RL4evTRRwPP2batYcOGqaamRm+//baee+45TZs2Tb/+9a9dv7+WKh39+v7776tDhw6aPn26Pv30U/3yl7/UhAkT9Mc//tH1+2up0tGvwW644Qb16dPHlXvBYenq1yeeeEK//OUvdc899+jTTz/VP//5Tw0ePNjVe2vp0tG3+/fv15AhQ9StWzetWLFCS5cuVVFRkQYPHqza2lrX77ElSkW/Ll++XEOGDNFFF12kd999VytXrtT48ePl8Rz+MXrUqFH69NNPtWDBAs2bN09vvfWWfvSjH6XlHlsMgyOSJPPKK68EHq9Zs8ZIMp988kngmG3bpn379uZ//ud/AsfOPfdcc9ttt8V839dee814PB6zdevWwLHJkyeb4uJiU11dndJ7QCS3+jWam266yQwcOLCxTUYS3O7XP//5z+bcc881CxcuNJLMnj17Uth6xOJWv+7evdvk5eWZf/7zn240G0lwq29XrlxpJJmNGzcGjn388cdGklm7dm1K7wGRGtqv/fv3N/fee2/M9/3ss8+MJLNy5crAsddff91YlmW+/vrr1N5EC8ZIS4aorq6WJOXm5gaOeTwe5eTkaOnSpSHnzpgxQ+3atdPJJ5+sCRMm6MCBA4Hnli9frlNOOUUdO3YMHBs8eLDKy8v16aefunwXCJeqfo1m3759atOmTeobjYRS2a+fffaZHnjgAf31r38N+a0f0i9V/bpgwQI5jqOvv/5avXv3VteuXTVixAht2rQpPTeCCKnq2169eqlt27Z65plnVFNTo4MHD+qZZ55R79691aNHj7TcCw5Lpl+3b9+uFStWqEOHDjrrrLPUsWNHnXvuuSH9vnz5cpWWlqpv376BY4MGDZLH49GKFSvSdDeZj//DZYgTTjhB3bp104QJE7Rnzx7V1NTod7/7nTZv3qwtW7YEzrv66qs1ffp0LVq0SBMmTNDf/vY3XXPNNYHnt27dGhJYJAUeb926NT03g4BU9Wu4t99+W7NmzWLouomkql+rq6s1cuRI/f73v1e3bt2a4lYQJFX9+u9//1uO4+i3v/2tnnrqKb300kvavXu3LrzwQtXU1DTFrbV4qerboqIiLV68WNOnT1deXp4KCwv1xhtv6PXXX1dWVlZT3FqLlky//vvf/5Yk3X///Ro3bpzeeOMNnX766brgggsCtS9bt25Vhw4dQt47KytLbdq04WenFOI7JEO0atVKf//733XDDTeoTZs28nq9GjRokIYOHSpjTOC84B9STznlFHXu3FkXXHCB1q1bp2OPPbYpmo443OjXTz75RMOHD9d9992niy66KG33gsNS1a8TJkxQ79694wZUpE+q+tVxHNXW1uoPf/hD4Hv0+eefV6dOnbRo0SJqW5pAqvr24MGDuuGGG/Sd73xHzz//vGzb1mOPPaZhw4Zp5cqVysvLa4rba7GS6VfHcSRJP/7xj3XddddJkk477TQtXLhQzz77rB5++OEma39Lw0hLBjnjjDO0atUq7d27V1u2bNEbb7yhXbt26Zhjjon5mv79+0uSvvrqK0lSp06dtG3btpBz/I87derkUssRTyr61e+zzz7TBRdcoB/96Ee69957XW034ktFv/7rX//Siy++qKysLGVlZemCCy6QJLVr10733Xef+zeBCKno186dO0uSTjzxxMA57du3V7t27bRx40YXW494UtG3M2fO1IYNGzR16lT169dP3/72tzVz5kytX79ec+fOTct9IFSifo32/ShJvXv3Dnw/durUSdu3bw953ufzaffu3fzslEKElgxUUlKi9u3ba+3atXrvvfc0fPjwmOeuWrVK0uFvygEDBmj16tUh33wLFixQcXFxxDcs0qsx/SpJn376qQYOHKgxY8booYcecru5SFJj+vXll1/WRx99pFWrVmnVqlV6+umnJUlLlixhedwm1ph+/c53viNJWrNmTeCc3bt3a+fOnerevbt7jUZSGtO3Bw4ckMfjkWVZgXP8j/2/0UfTiNWvPXr0UJcuXUK+HyXpyy+/DHw/DhgwQHv37tX7778feP5f//qXHMcJBFekQFOuAoD6qaioMB9++KH58MMPjSTzxBNPmA8//ND85z//McYYM3v2bLNo0SKzbt06M2fOHNO9e3dzxRVXBF7/1VdfmQceeMC89957Zv369Wbu3LnmmGOOMeecc07gHJ/PZ04++WRz0UUXmVWrVpk33njDtG/f3kyYMCHt99tSpKNfV69ebdq3b2+uueYas2XLlsDX9u3b036/LUU6+jXcokWLWD3MZenq1+HDh5uTTjrJLFu2zKxevdpccskl5sQTTzQ1NTVpvd+WJB19+/nnn5ucnBxz4403ms8++8x88skn5pprrjElJSXmm2++Sfs9twSN7VdjjHnyySdNcXGxefHFF83atWvNvffea3Jzc81XX30VOGfIkCHmtNNOMytWrDBLly41xx13nBk5cmRa7zXTEVqOIP4fSMK/xowZY4wxZtKkSaZr166mVatWplu3bubee+8NWaZ448aN5pxzzjFt2rQxOTk5pmfPnuauu+4y+/btC7nOhg0bzNChQ01eXp5p166d+dnPfmZqa2vTeastSjr69b777ot6je7du6f5bluOdH2/RrsmocU96erXffv2meuvv96UlpaaNm3amMsvvzxkmVykXrr69v/+7//Md77zHVNSUmJat25tzj//fLN8+fJ03mqL0th+9Xv44YdN165dTX5+vhkwYIBZsmRJyPO7du0yI0eONIWFhaa4uNhcd911pqKiIh232GJYxgRVkAEAAABAM0NNCwAAAIBmjdACAAAAoFkjtAAAAABo1ggtAAAAAJo1QgsAAACAZo3QAgAAAKBZI7QAAAAAaNYILQAAAACaNUILAAAAgGaN0AIAAACgWSO0AAAAAGjWCC0AAAAAmrX/D/REKaJdW53NAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import calendar\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from river import compose\n",
    "from river import datasets\n",
    "from river import linear_model\n",
    "from river import metrics\n",
    "from river import optim\n",
    "from river import preprocessing\n",
    "from river import stats\n",
    "    \n",
    "\n",
    "def get_ordinal_date(x):\n",
    "    return {'ordinal_date': x['month'].toordinal()}    \n",
    "\n",
    "    \n",
    "def get_month_distances(x):\n",
    "    return {\n",
    "        calendar.month_name[month]: math.exp(-(x['month'].month - month) ** 2)\n",
    "        for month in range(1, 13)\n",
    "    }\n",
    "        \n",
    "\n",
    "def make_model(alpha):\n",
    "    \n",
    "    extract_features = compose.TransformerUnion(get_ordinal_date, get_month_distances)\n",
    "\n",
    "    scale = preprocessing.StandardScaler()\n",
    "\n",
    "    learn = linear_model.LinearRegression(\n",
    "        intercept_lr=0,\n",
    "        optimizer=optim.SGD(0.03),\n",
    "        loss=optim.losses.Quantile(alpha=alpha)\n",
    "    )\n",
    "\n",
    "    model = extract_features | scale | learn\n",
    "    model = preprocessing.TargetStandardScaler(regressor=model)\n",
    "\n",
    "    return model\n",
    "\n",
    "metric = metrics.MAE()\n",
    "\n",
    "models = {\n",
    "    'lower': make_model(alpha=0.05),\n",
    "    'center': make_model(alpha=0.5),\n",
    "    'upper': make_model(alpha=0.95)\n",
    "}\n",
    "\n",
    "dates = []\n",
    "y_trues = []\n",
    "y_preds = {\n",
    "    'lower': [],\n",
    "    'center': [],\n",
    "    'upper': []\n",
    "}\n",
    "\n",
    "for x, y in datasets.AirlinePassengers():\n",
    "    y_trues.append(y)\n",
    "    dates.append(x['month'])\n",
    "    \n",
    "    for name, model in models.items():\n",
    "        y_preds[name].append(model.predict_one(x))\n",
    "        model.learn_one(x, y)\n",
    "\n",
    "    # Update the error metric\n",
    "    metric.update(y, y_preds['center'][-1])\n",
    "\n",
    "# Plot the results\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.grid(alpha=0.75)\n",
    "ax.plot(dates, y_trues, lw=3, color='#2ecc71', alpha=0.8, label='Truth')\n",
    "ax.plot(dates, y_preds['center'], lw=3, color='#e74c3c', alpha=0.8, label='Prediction')\n",
    "ax.fill_between(dates, y_preds['lower'], y_preds['upper'], color='#e74c3c', alpha=0.3, label='Prediction interval')\n",
    "ax.legend()\n",
    "ax.set_title(metric);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important thing to note is that the prediction interval we obtained should not be confused with a confidence interval. Simply put, a prediction interval represents uncertainty for where the true value lies, whereas a confidence interval encapsulates the uncertainty on the prediction. You can find out more by reading [this](https://stats.stackexchange.com/questions/16493/difference-between-confidence-intervals-and-prediction-intervals) CrossValidated post."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
